#define fast_norm3( v ) (v).w = native_rsqrt((v).x*(v).x + (v).y*(v).y + (v).z+(v).z); (v).x *= (v).w; (v).y *= (v).w; (v).z *= (v).w; (poot).w = 0.0f;
#define fast_norm3_t( temp, v ) (temp) = native_rsqrt((v).x*(v).x + (v).y*(v).y + (v).z*(v).z ); (v).x *= (temp); (v).y *= (temp); (v).z *= (temp);

#define fast_dot3( v1, v2 ) ((v1).x*(v2).x + (v1).y*(v2).y + (v1).z*(v2).z)
#define fast_sub3( r, v1, v2 ) (r).x = (v1).x - (v2).x; (r).y = (v1).y - (v2).y; (r).z = (v1).z - (v2).z;
#define fast_add3( r, v1, v2 ) (r).x = (v1).x + (v2).x; (r).y = (v1).y + (v2).y; (r).z = (v1).z + (v2).z;
#define fast_min3( r, v1, v2 ) (r).x = fmin((v1).x, (v2).x); (r).y = fmin((v1).y, (v2).y); (r).z = fmin((v1).z, (v2).z);
#define fast_max3( r, v1, v2 ) (r).x = fmax((v1).x, (v2).x); (r).y = fmax((v1).y, (v2).y); (r).z = fmax((v1).z, (v2).z);

typedef struct {
	float x;
	float y;
	float z;
} float3h;

#define NUM_THREADS ((BLOCK_X) * (BLOCK_Y))

#define REDUCE_MIN(i) \
if ( threadID < (i) ) { \
	value = fmin(value, array[threadID + (i)]); \
	array[threadID] = value;\
}\
barrier(CLK_LOCAL_MEM_FENCE);

#define REDUCE_MAX(i) \
if ( threadID < (i) ) { \
	value = fmax(value, array[threadID + (i)]); \
	array[threadID] = value;\
}\
barrier(CLK_LOCAL_MEM_FENCE);

#define REDUCE_ALL(i) \
if ( threadID < (i) ) { \
	normal.x += array[(i) + (0 * NUM_THREADS) ]; \
	normal.y += array[(i) + (1 * NUM_THREADS) ]; \
	normal.z += array[(i) + (2 * NUM_THREADS) ]; \
	minimum = fmin(minimum, array[(i) + (3 * NUM_THREADS)]); \
	maximum = fmax(maximum, array[(i) + (4 * NUM_THREADS)]); \
	array[(0 * NUM_THREADS)] = normal.x; \
	array[(1 * NUM_THREADS)] = normal.y; \
	array[(2 * NUM_THREADS)] = normal.z; \
	array[(3 * NUM_THREADS)] = minimum; \
	array[(4 * NUM_THREADS)] = maximum; \
}\
barrier(CLK_LOCAL_MEM_FENCE);

#define REDUCE_AVERAGE(i) \
if ( threadID < (i) ) { \
	value.x += array[threadID + (i) + 0 * NUM_THREADS ]; \
	value.y += array[threadID + (i) + 1 * NUM_THREADS ]; \
	value.z += array[threadID + (i) + 2 * NUM_THREADS ]; \
	array[(0 * NUM_THREADS) + threadID] = value.x; \
	array[(1 * NUM_THREADS) + threadID] = value.y; \
	array[(2 * NUM_THREADS) + threadID] = value.z; \
}\
barrier(CLK_LOCAL_MEM_FENCE);


void reduceAll( float3h normal, float minimum, float maximum, __local float *array, int threadID) {

	array[0 * NUM_THREADS] = normal.x;
	array[1 * NUM_THREADS] = normal.y;
	array[2 * NUM_THREADS] = normal.z;
	array[3 * NUM_THREADS] = minimum;
	array[4 * NUM_THREADS] = maximum;
	barrier(CLK_LOCAL_MEM_FENCE);
	
	#if NUM_THREADS > 128
	REDUCE_ALL(128);
	#endif
	#if NUM_THREADS > 64
	REDUCE_ALL(64);
	#endif
	#if NUM_THREADS > 32
	REDUCE_ALL(32);
	#endif
	REDUCE_ALL(16);
	REDUCE_ALL(8);
	REDUCE_ALL(4);
	REDUCE_ALL(2);
	REDUCE_ALL(1);
	
}


float3h reduceAdd( float3h value, __local float *array, int threadID) {

	array[0 * NUM_THREADS + threadID] = value.x;
	array[1 * NUM_THREADS + threadID] = value.y;
	array[2 * NUM_THREADS + threadID] = value.z;
	barrier(CLK_LOCAL_MEM_FENCE);
	
	#if NUM_THREADS > 128
	REDUCE_AVERAGE(128);
	#endif
	#if NUM_THREADS > 64
	REDUCE_AVERAGE(64);
	#endif
	#if NUM_THREADS > 32
	REDUCE_AVERAGE(32);
	#endif
	REDUCE_AVERAGE(16);
	REDUCE_AVERAGE(8);
	REDUCE_AVERAGE(4);
	REDUCE_AVERAGE(2);
	REDUCE_AVERAGE(1);

	float3h rValue;
	rValue.x = array[0 * NUM_THREADS];
	rValue.y = array[1 * NUM_THREADS];
	rValue.z = array[2 * NUM_THREADS];
	return rValue;
	
}

float reduceMax( float value, __local float *array, int threadID) {

	array[threadID] = value;
	barrier(CLK_LOCAL_MEM_FENCE);

	#if NUM_THREADS > 128
	REDUCE_MAX(128)
	#endif
	#if NUM_THREADS > 64
	REDUCE_MAX(64)
	#endif
	#if NUM_THREADS > 32
	REDUCE_MAX(32)
	#endif
	REDUCE_MAX(16)
	REDUCE_MAX(8)
	REDUCE_MAX(4)
	REDUCE_MAX(2)
	REDUCE_MAX(1)

	return array[0];

}


float reduceMin( float value, __local float *array, int threadID) {

	array[threadID] = value;
	barrier(CLK_LOCAL_MEM_FENCE);

	#if NUM_THREADS > 128
	REDUCE_MIN(128);
	#endif
	#if NUM_THREADS > 64
	REDUCE_MIN(64);
	#endif
	#if NUM_THREADS > 32
	REDUCE_MIN(32);
	#endif
	REDUCE_MIN(16);
	REDUCE_MIN(8);
	REDUCE_MIN(4);
	REDUCE_MIN(2);
	REDUCE_MIN(1);

	return array[0];

}

__kernel void perblockattributes(	__read_only	 image2d_t	normSpecImage,
									__read_only	 image2d_t	depthImage,
									const		sampler_t	nearestSampler,
									const		int2		imageSize,
									const		int2		num_blocks,
									__write_only image2d_t	perBlockImage1, 
									__write_only image2d_t	perBlockImage2)
{
	
	float temp;
	__local float reduceArray[ NUM_THREADS * 5 ];


	const int2 gid = {	get_global_id(0), get_global_id(1) };
						
	const int linearThreadID = get_local_id(0) + get_local_size(0) * get_local_id(1);

	float3h normal;
	float depth;
	
	bool validPixel = (gid.x < imageSize.x && gid.y < imageSize.y);
	
	if ( validPixel ) {
	
		// get the normal, normalize it
		float4	normal_lookup	= read_imagef(normSpecImage, nearestSampler, gid);
		normal.x = normal_lookup.x - 0.5f;
		normal.y = normal_lookup.y - 0.5f;
		normal.z = normal_lookup.z - 0.5f;
		fast_norm3_t(temp, normal);

		// get the depth
		depth = read_imagef(depthImage, nearestSampler, gid).x;
		
	}
	else {
		// don't affect sum of normals
		normal.x = 0.0f;
		normal.y = 0.0f;
		normal.z = 0.0f;
	}
	
	
	// if it's not a valid pixel, set up the min/max depth so as to
	// not affect block computation
	float pMinZ			= validPixel ? depth : 1.0f;
	float pMaxZ			= validPixel ? depth : 0.0f;
	
	reduceAll( normal, pMinZ, pMaxZ, &(reduceArray[linearThreadID]), linearThreadID );
	float minDepth		= reduceArray[( 3 * NUM_THREADS )];	
	float maxDepth		= reduceArray[( 4 * NUM_THREADS )];;
	
	float3h	average_normal;
	average_normal.x		= reduceArray[( 0 * NUM_THREADS )];	
	average_normal.y		= reduceArray[( 1 * NUM_THREADS )];
	average_normal.z		= reduceArray[( 2 * NUM_THREADS )];
	fast_norm3_t(temp, average_normal);
	
	// don't affect normal cone size if not valid pixel
	float normalDot		= validPixel ? fast_dot3(normal, average_normal) : 1.0f;
	float minDot		= reduceMin( minDot, reduceArray, linearThreadID);
	
	if ( linearThreadID == 0 ) {
	
		int2 bid = { get_group_id(0), get_group_id(1) };
	
		float normal_cone_size = minDot <= 0.0f ? 100.0f : 1.0f - minDot * minDot;

		/* compute ndc versions of min and max depths */
		float maxZ = ( 2.0f * maxDepth - 1.0f );
		float minZ = ( 2.0f * minDepth - 1.0f );	
		
		normal_cone_size = 100.0f;
								
		float4 out1 = { average_normal.x, average_normal.y, average_normal.z, normal_cone_size };
		float4 out2 = { minZ, maxZ, 0.0f, 0.0f };

		write_imagef( perBlockImage1, bid, out1 );
		write_imagef( perBlockImage2, bid, out2 );
		
	}

}